<!DOCTYPE html>

<html xmlns="http://www.w3.org/1999/xhtml">

<head>

<meta charset="utf-8" />
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="pandoc" />

<meta name="viewport" content="width=device-width, initial-scale=1">

<meta name="author" content="Florian P Breitwieser" />

<meta name="date" content="2017-11-05" />

<title>Using pavian in R</title>



<style type="text/css">code{white-space: pre;}</style>
<style type="text/css">
div.sourceCode { overflow-x: auto; }
table.sourceCode, tr.sourceCode, td.lineNumbers, td.sourceCode {
  margin: 0; padding: 0; vertical-align: baseline; border: none; }
table.sourceCode { width: 100%; line-height: 100%; }
td.lineNumbers { text-align: right; padding-right: 4px; padding-left: 4px; color: #aaaaaa; border-right: 1px solid #aaaaaa; }
td.sourceCode { padding-left: 5px; }
code > span.kw { color: #007020; font-weight: bold; } /* Keyword */
code > span.dt { color: #902000; } /* DataType */
code > span.dv { color: #40a070; } /* DecVal */
code > span.bn { color: #40a070; } /* BaseN */
code > span.fl { color: #40a070; } /* Float */
code > span.ch { color: #4070a0; } /* Char */
code > span.st { color: #4070a0; } /* String */
code > span.co { color: #60a0b0; font-style: italic; } /* Comment */
code > span.ot { color: #007020; } /* Other */
code > span.al { color: #ff0000; font-weight: bold; } /* Alert */
code > span.fu { color: #06287e; } /* Function */
code > span.er { color: #ff0000; font-weight: bold; } /* Error */
code > span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
code > span.cn { color: #880000; } /* Constant */
code > span.sc { color: #4070a0; } /* SpecialChar */
code > span.vs { color: #4070a0; } /* VerbatimString */
code > span.ss { color: #bb6688; } /* SpecialString */
code > span.im { } /* Import */
code > span.va { color: #19177c; } /* Variable */
code > span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code > span.op { color: #666666; } /* Operator */
code > span.bu { } /* BuiltIn */
code > span.ex { } /* Extension */
code > span.pp { color: #bc7a00; } /* Preprocessor */
code > span.at { color: #7d9029; } /* Attribute */
code > span.do { color: #ba2121; font-style: italic; } /* Documentation */
code > span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code > span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code > span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
</style>



<link href="data:text/css;charset=utf-8,body%20%7B%0Abackground%2Dcolor%3A%20%23fff%3B%0Amargin%3A%201em%20auto%3B%0Amax%2Dwidth%3A%20700px%3B%0Aoverflow%3A%20visible%3B%0Apadding%2Dleft%3A%202em%3B%0Apadding%2Dright%3A%202em%3B%0Afont%2Dfamily%3A%20%22Open%20Sans%22%2C%20%22Helvetica%20Neue%22%2C%20Helvetica%2C%20Arial%2C%20sans%2Dserif%3B%0Afont%2Dsize%3A%2014px%3B%0Aline%2Dheight%3A%201%2E35%3B%0A%7D%0A%23header%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0A%23TOC%20%7B%0Aclear%3A%20both%3B%0Amargin%3A%200%200%2010px%2010px%3B%0Apadding%3A%204px%3B%0Awidth%3A%20400px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Aborder%2Dradius%3A%205px%3B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Afont%2Dsize%3A%2013px%3B%0Aline%2Dheight%3A%201%2E3%3B%0A%7D%0A%23TOC%20%2Etoctitle%20%7B%0Afont%2Dweight%3A%20bold%3B%0Afont%2Dsize%3A%2015px%3B%0Amargin%2Dleft%3A%205px%3B%0A%7D%0A%23TOC%20ul%20%7B%0Apadding%2Dleft%3A%2040px%3B%0Amargin%2Dleft%3A%20%2D1%2E5em%3B%0Amargin%2Dtop%3A%205px%3B%0Amargin%2Dbottom%3A%205px%3B%0A%7D%0A%23TOC%20ul%20ul%20%7B%0Amargin%2Dleft%3A%20%2D2em%3B%0A%7D%0A%23TOC%20li%20%7B%0Aline%2Dheight%3A%2016px%3B%0A%7D%0Atable%20%7B%0Amargin%3A%201em%20auto%3B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dcolor%3A%20%23DDDDDD%3B%0Aborder%2Dstyle%3A%20outset%3B%0Aborder%2Dcollapse%3A%20collapse%3B%0A%7D%0Atable%20th%20%7B%0Aborder%2Dwidth%3A%202px%3B%0Apadding%3A%205px%3B%0Aborder%2Dstyle%3A%20inset%3B%0A%7D%0Atable%20td%20%7B%0Aborder%2Dwidth%3A%201px%3B%0Aborder%2Dstyle%3A%20inset%3B%0Aline%2Dheight%3A%2018px%3B%0Apadding%3A%205px%205px%3B%0A%7D%0Atable%2C%20table%20th%2C%20table%20td%20%7B%0Aborder%2Dleft%2Dstyle%3A%20none%3B%0Aborder%2Dright%2Dstyle%3A%20none%3B%0A%7D%0Atable%20thead%2C%20table%20tr%2Eeven%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Ap%20%7B%0Amargin%3A%200%2E5em%200%3B%0A%7D%0Ablockquote%20%7B%0Abackground%2Dcolor%3A%20%23f6f6f6%3B%0Apadding%3A%200%2E25em%200%2E75em%3B%0A%7D%0Ahr%20%7B%0Aborder%2Dstyle%3A%20solid%3B%0Aborder%3A%20none%3B%0Aborder%2Dtop%3A%201px%20solid%20%23777%3B%0Amargin%3A%2028px%200%3B%0A%7D%0Adl%20%7B%0Amargin%2Dleft%3A%200%3B%0A%7D%0Adl%20dd%20%7B%0Amargin%2Dbottom%3A%2013px%3B%0Amargin%2Dleft%3A%2013px%3B%0A%7D%0Adl%20dt%20%7B%0Afont%2Dweight%3A%20bold%3B%0A%7D%0Aul%20%7B%0Amargin%2Dtop%3A%200%3B%0A%7D%0Aul%20li%20%7B%0Alist%2Dstyle%3A%20circle%20outside%3B%0A%7D%0Aul%20ul%20%7B%0Amargin%2Dbottom%3A%200%3B%0A%7D%0Apre%2C%20code%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0Aborder%2Dradius%3A%203px%3B%0Acolor%3A%20%23333%3B%0Awhite%2Dspace%3A%20pre%2Dwrap%3B%20%0A%7D%0Apre%20%7B%0Aborder%2Dradius%3A%203px%3B%0Amargin%3A%205px%200px%2010px%200px%3B%0Apadding%3A%2010px%3B%0A%7D%0Apre%3Anot%28%5Bclass%5D%29%20%7B%0Abackground%2Dcolor%3A%20%23f7f7f7%3B%0A%7D%0Acode%20%7B%0Afont%2Dfamily%3A%20Consolas%2C%20Monaco%2C%20%27Courier%20New%27%2C%20monospace%3B%0Afont%2Dsize%3A%2085%25%3B%0A%7D%0Ap%20%3E%20code%2C%20li%20%3E%20code%20%7B%0Apadding%3A%202px%200px%3B%0A%7D%0Adiv%2Efigure%20%7B%0Atext%2Dalign%3A%20center%3B%0A%7D%0Aimg%20%7B%0Abackground%2Dcolor%3A%20%23FFFFFF%3B%0Apadding%3A%202px%3B%0Aborder%3A%201px%20solid%20%23DDDDDD%3B%0Aborder%2Dradius%3A%203px%3B%0Aborder%3A%201px%20solid%20%23CCCCCC%3B%0Amargin%3A%200%205px%3B%0A%7D%0Ah1%20%7B%0Amargin%2Dtop%3A%200%3B%0Afont%2Dsize%3A%2035px%3B%0Aline%2Dheight%3A%2040px%3B%0A%7D%0Ah2%20%7B%0Aborder%2Dbottom%3A%204px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Apadding%2Dbottom%3A%202px%3B%0Afont%2Dsize%3A%20145%25%3B%0A%7D%0Ah3%20%7B%0Aborder%2Dbottom%3A%202px%20solid%20%23f7f7f7%3B%0Apadding%2Dtop%3A%2010px%3B%0Afont%2Dsize%3A%20120%25%3B%0A%7D%0Ah4%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23f7f7f7%3B%0Amargin%2Dleft%3A%208px%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Ah5%2C%20h6%20%7B%0Aborder%2Dbottom%3A%201px%20solid%20%23ccc%3B%0Afont%2Dsize%3A%20105%25%3B%0A%7D%0Aa%20%7B%0Acolor%3A%20%230033dd%3B%0Atext%2Ddecoration%3A%20none%3B%0A%7D%0Aa%3Ahover%20%7B%0Acolor%3A%20%236666ff%3B%20%7D%0Aa%3Avisited%20%7B%0Acolor%3A%20%23800080%3B%20%7D%0Aa%3Avisited%3Ahover%20%7B%0Acolor%3A%20%23BB00BB%3B%20%7D%0Aa%5Bhref%5E%3D%22http%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0Aa%5Bhref%5E%3D%22https%3A%22%5D%20%7B%0Atext%2Ddecoration%3A%20underline%3B%20%7D%0A%0Acode%20%3E%20span%2Ekw%20%7B%20color%3A%20%23555%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Edt%20%7B%20color%3A%20%23902000%3B%20%7D%20%0Acode%20%3E%20span%2Edv%20%7B%20color%3A%20%2340a070%3B%20%7D%20%0Acode%20%3E%20span%2Ebn%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Efl%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Ech%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Est%20%7B%20color%3A%20%23d14%3B%20%7D%20%0Acode%20%3E%20span%2Eco%20%7B%20color%3A%20%23888888%3B%20font%2Dstyle%3A%20italic%3B%20%7D%20%0Acode%20%3E%20span%2Eot%20%7B%20color%3A%20%23007020%3B%20%7D%20%0Acode%20%3E%20span%2Eal%20%7B%20color%3A%20%23ff0000%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%0Acode%20%3E%20span%2Efu%20%7B%20color%3A%20%23900%3B%20font%2Dweight%3A%20bold%3B%20%7D%20%20code%20%3E%20span%2Eer%20%7B%20color%3A%20%23a61717%3B%20background%2Dcolor%3A%20%23e3d2d2%3B%20%7D%20%0A" rel="stylesheet" type="text/css" />

</head>

<body>




<h1 class="title toc-ignore">Using pavian in R</h1>
<h4 class="author"><em>Florian P Breitwieser</em></h4>
<h4 class="date"><em>2017-11-05</em></h4>



<p>Most functions that are available in the Pavian shiny interface can also be used from the command line.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">library</span>(pavian)
sample_data &lt;-<span class="st"> </span>pavian::<span class="kw">read_sample_data</span>(<span class="kw">system.file</span>(<span class="st">&quot;shinyapp&quot;</span>,<span class="st">&quot;example-data&quot;</span>,<span class="st">&quot;brain-biopsies&quot;</span>,<span class="dt">package=</span><span class="st">&quot;pavian&quot;</span>))
reports &lt;-<span class="st"> </span>pavian::<span class="kw">read_reports</span>(sample_data$ReportFilePath, sample_data$Name)</code></pre></div>
<pre><code>## [Nov05 10:44] Reading 11 reports ...</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">merged_reports &lt;-<span class="st"> </span>pavian::<span class="kw">merge_reports2</span>(reports, <span class="dt">col_names =</span> sample_data$Name)

pavian::<span class="kw">summarize_reports</span>(reports)</code></pre></div>
<pre><code>##        number_of_raw_reads classified_reads chordate_reads
## PT1               12022284         11951847       11245467
## PT2                8294101          8281974        8219504
## PT3               17669644         17609799       16992534
## PT4               29101779         29038413       28453071
## PT5               26919065         26715453       25903121
## PT6               27261739         27159918       26654128
## PT7               19065574         18918495       13462479
## PT8-S1            14957133         14925304       14643350
## PT8-S2            13990253         13961725       13952205
## PT9               26500914         26422494       26183265
## PT10              21319274         21244417       21027311
##        artificial_reads unclassified_reads microbial_reads bacterial_reads
## PT1              704383              70437            1790            1497
## PT2               62140              12127             304             141
## PT3              610050              59845            7099            6859
## PT4              583338              63366            1821            1094
## PT5              786873             203612           25314           16046
## PT6              502575             101821            3034            2854
## PT7             5446842             147079            9089            8781
## PT8-S1           281657              31829              58              18
## PT8-S2             9225              28528              52              22
## PT9              236228              78420            2695            2052
## PT10             215939              74857            1029             569
##        viral_reads fungal_reads protozoan_reads
## PT1              0           73              27
## PT2              0            0               0
## PT3              3           58               4
## PT4              0           15              21
## PT5           8954           50               8
## PT6             17           11               3
## PT7             74           13               4
## PT8-S1           1            6               1
## PT8-S2           1            1               0
## PT9              1          519               0
## PT10            20          180               0</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">tax_data &lt;-<span class="st"> </span>merged_reports[[<span class="dv">1</span>]]
clade_reads &lt;-<span class="st"> </span>merged_reports[[<span class="dv">2</span>]]
taxon_reads &lt;-<span class="st"> </span>merged_reports[[<span class="dv">3</span>]]

<span class="kw">colSums</span>(clade_reads,<span class="dt">na.rm =</span> T)</code></pre></div>
<pre><code>##       PT1       PT2       PT3       PT4       PT5       PT6       PT7 
## 362757256 263286609 546326888 912906520 832441890 855069673 452821122 
##    PT8-S1    PT8-S2       PT9      PT10 
## 469746275 446536577 838915270 673820947</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw">colSums</span>(taxon_reads,<span class="dt">na.rm =</span> T)</code></pre></div>
<pre><code>##      PT1      PT2      PT3      PT4      PT5      PT6      PT7   PT8-S1 
## 12022284  8294101 17669644 29101779 26919065 27261739 19065574 14957133 
##   PT8-S2      PT9     PT10 
## 13990253 26500914 21319274</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">sel_rows &lt;-<span class="st"> </span>pavian::<span class="kw">filter_taxa</span>(tax_data,
                                <span class="dt">rm_clades =</span> <span class="kw">c</span>(<span class="st">&quot;Chordata&quot;</span>, <span class="st">&quot;artificial sequences&quot;</span>, <span class="st">&quot;unclassified&quot;</span>),
                                <span class="dt">taxRank =</span> <span class="st">&quot;S&quot;</span>)
<span class="kw">summary</span>(sel_rows)</code></pre></div>
<pre><code>##    Mode   FALSE    TRUE 
## logical    2339    1043</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">filtered_clade_reads &lt;-<span class="st"> </span>pavian::<span class="kw">filter_cladeReads</span>(clade_reads, tax_data, <span class="kw">c</span>(<span class="st">&quot;Chordata&quot;</span>, <span class="st">&quot;artificial sequences&quot;</span>, <span class="st">&quot;unclassified&quot;</span>))

tax_data1 &lt;-<span class="st"> </span>tax_data[sel_rows,]
filtered_clade_reads1 &lt;-<span class="st"> </span>filtered_clade_reads[sel_rows, ]
taxon_reads1 &lt;-<span class="st"> </span>taxon_reads[sel_rows, ]

<span class="kw">head</span>(<span class="kw">cbind</span>(tax_data1[,<span class="dv">1</span>:<span class="dv">3</span>],clade_reads[sel_rows, ])[<span class="kw">order</span>(-<span class="kw">apply</span>(filtered_clade_reads1,<span class="dv">1</span>,max, <span class="dt">na.rm=</span>T)),])</code></pre></div>
<pre><code>##                                name taxRank  taxID PT1 PT2 PT3 PT4  PT5
## 2804                JC polyomavirus       S  10632  NA  NA  NA  NA 8944
## 406         Propionibacterium acnes       S   1747   2  15  61  87  539
## 1937        Akkermansia muciniphila       S 239935  NA  NA 121   1 1031
## 113                Escherichia coli       S    562 750  44 640 250  548
## 2148 Candidatus Nitrospira defluvii       S 330214  NA  NA   7  NA  712
## 2354 Candidatus Babela massiliensis       S 673862  NA  NA  NA  NA  659
##      PT6  PT7 PT8-S1 PT8-S2 PT9 PT10
## 2804  15   NA     NA     NA  NA   NA
## 406   63 1459      2      3 222   23
## 1937 102    7     NA     NA  14   NA
## 113  447  148     NA     NA 174  190
## 2148  12    6     NA     NA  NA    1
## 2354   7    1     NA     NA  NA   NA</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">normalized_clade_reads &lt;-<span class="st"> </span><span class="kw">normalize</span>(filtered_clade_reads1)
normalized_taxon_reads &lt;-<span class="st"> </span><span class="kw">normalize</span>(taxon_reads[sel_rows,], <span class="dt">sums =</span> <span class="kw">colSums</span>(filtered_clade_reads1,<span class="dt">na.rm =</span> T))
<span class="kw">head</span>(<span class="kw">cbind</span>(tax_data1[,<span class="dv">1</span>:<span class="dv">3</span>],<span class="dt">max=</span><span class="kw">apply</span>(<span class="kw">cbind</span>(normalized_clade_reads),<span class="dv">1</span>,max, <span class="dt">na.rm=</span>T), normalized_clade_reads)[<span class="kw">order</span>(-<span class="kw">apply</span>(<span class="kw">cbind</span>(normalized_clade_reads),<span class="dv">1</span>,max, <span class="dt">na.rm=</span>T)),])</code></pre></div>
<pre><code>##                          name taxRank taxID       max         PT1
## 113          Escherichia coli       S   562 0.8134490 0.813449024
## 2804          JC polyomavirus       S 10632 0.4929997          NA
## 55   Saccharomyces cerevisiae       S  4932 0.3392857 0.003253796
## 406   Propionibacterium acnes       S  1747 0.3000000 0.002169197
## 213       Delftia acidovorans       S 80866 0.1363636 0.003253796
## 172        Pseudomonas putida       S   303 0.1000000 0.004338395
##             PT2          PT3        PT4          PT5         PT6
## 113  0.43564356 0.1236954001 0.43554007 3.020615e-02 0.295634921
## 2804         NA           NA         NA 4.929997e-01 0.009920635
## 55           NA 0.0001932741 0.00174216 5.512071e-05          NA
## 406  0.14851485 0.0117897178 0.15156794 2.971007e-02 0.041666667
## 213  0.08910891 0.0168148434         NA           NA          NA
## 172  0.02970297 0.0094704291         NA 2.204829e-04 0.005291005
##              PT7     PT8-S1 PT8-S2         PT9       PT10
## 113  0.020029774         NA     NA 0.086739781 0.37698413
## 2804          NA         NA     NA          NA         NA
## 55            NA 0.22727273    0.1 0.251744766 0.33928571
## 406  0.197455677 0.09090909    0.3 0.110667996 0.04563492
## 213  0.017323048 0.13636364     NA 0.001994018         NA
## 172  0.008526188 0.04545455    0.1 0.004486540         NA</code></pre>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r">reads_zscore &lt;-<span class="st"> </span><span class="kw">robust_zscore</span>(<span class="dv">100</span>*<span class="kw">cbind</span>(normalized_clade_reads,normalized_taxon_reads), <span class="fl">0.001</span>)
clade_reads_zscore &lt;-<span class="st"> </span>reads_zscore[,<span class="dv">1</span>:<span class="dv">20</span>]
reads_zscore_df &lt;-<span class="st"> </span><span class="kw">cbind</span>(tax_data1[,<span class="dv">1</span>:<span class="dv">3</span>],<span class="dt">max=</span><span class="kw">apply</span>(clade_reads_zscore,<span class="dv">1</span>,max, <span class="dt">na.rm=</span>T), clade_reads_zscore)[<span class="kw">order</span>(-<span class="kw">apply</span>(clade_reads_zscore,<span class="dv">1</span>,max, <span class="dt">na.rm=</span>T)),]
## Calculate z-score from the clade reads</code></pre></div>



<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    script.src  = "https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML";
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>

</body>
</html>
